set matsize 10000
set more off
set seed 22242

global path "****"
global output "****"

cap log close
log using "$output/power_analysis_102619.log", replace


foreach complier_scalar of numlist 2(1)4{

foreach percent_effect of numlist 0.01(0.02)0.07{

foreach age of numlist 20(5)60{

di "Age: " `age'
di "Complier Scaler: " `complier_scalar'
di "Percent Effect: " `percent_effect'


use "$path/mort_pop_by_age2016_nofull.dta", clear

gen dead16_2= dead16*2
qui summ dead16_2 [fw=tin] if age2016 >= `age' & age2016 < 65	

local y_all r(mean)

if `age' == 20{
local beta 0.2487829
local N 4152374
}
else if `age' == 25{
local beta 0.251674
local N 3698553
}
else if `age' == 30{
local beta 0.2653685 
local N 2994294
}
else if `age' == 35{
local beta 0.2947112
local N 2321377
}
else if `age' == 40{
local beta 0.3419995
local N 1791670
}
else if `age' == 45{
local beta 0.3575606    
local N 1358983
}
else if `age' == 50{
local beta 0.371675 
local N 937776
}
else if `age' == 55{
local beta .4311927
local N 551814
}
else if `age' == 60{
local beta  0.4883639 
local N 221452
}

local y_comp `complier_scalar'*`y_all'
local N_t  0.86*`N'
local N_c  0.14*`N'

local delta  `percent_effect'

local gamma `delta'*`y_comp'

global p_c = `y_all'
global p_t = `y_all' - `gamma'*`beta'

clear

qui set obs `N'

qui gen treatment = 1 if _n <= `N_t'

qui replace treatment = 0 if treatment != 1

if `percent_effect' == 0.01{
local perc_mat_name 01
}
else if `percent_effect' == 0.03{
local perc_mat_name 03
}
else if `percent_effect' == 0.05{
local perc_mat_name 05
}
else if `percent_effect' == 0.07{
local perc_mat_name 07
}


cap mat drop p_val_`complier_scalar'_`perc_mat_name'_`age'
foreach i of numlist 1/1000{
cap drop dead
qui gen dead = rbinomial(1,${p_t}) if treatment == 1
qui replace dead = rbinomial(1,${p_c}) if treatment == 0
qui ttest dead, by(treatment)
matrix p_val_`complier_scalar'_`perc_mat_name'_`age' = nullmat(p_val_`complier_scalar'_`perc_mat_name'_`age') \ r(p) 
}
}
}
}

drop treatment dead
set obs 1000

foreach complier_scalar of numlist 2(1)4{
foreach percent_effect of numlist 0.01(0.02)0.07{
foreach age of numlist 20(5)60{
if `percent_effect' == 0.01{
local perc_mat_name 01
}
else if `percent_effect' == 0.03{
local perc_mat_name 03
}
else if `percent_effect' == 0.05{
local perc_mat_name 05
}
else if `percent_effect' == 0.07{
local perc_mat_name 07
}
svmat p_val_`complier_scalar'_`perc_mat_name'_`age'
gen reject_`complier_scalar'_`perc_mat_name'_`age' = (p_val_`complier_scalar'_`perc_mat_name'_`age' < 0.05)
}
}
}


save "$output/power_simulation_results_102619.dta", replace


use "$output/power_simulation_results_102619.dta", clear

collapse (mean) reject* if p_val_2_01_201 != .
gen id = 1

reshape long reject_2_01_ reject_2_03_  reject_2_05_  reject_2_07_  reject_3_01_  reject_3_03_  reject_3_05_  reject_3_07_  ///
 reject_4_01_ reject_4_03_  reject_4_05_  reject_4_07_, i(id) j(reject)

 drop id

gen id = _n

reshape long reject_2_01_ reject_2_03_  reject_2_05_  reject_2_07_  reject_3_01_  reject_3_03_  reject_3_05_  reject_3_07_  ///
 reject_4_01_ reject_4_03_  reject_4_05_  reject_4_07_, i(id) j(reject)
 

rename reject age

foreach i in 2_01 2_03 2_05 2_07 3_01 3_03 3_05 3_07 4_01 4_03 4_05 4_07{
	gen power`i' = reject_`i'_*100
	drop reject_`i'_
}  

drop id

save "$output/power_results_by_age_102619.dta", replace
